February 1, 2008 



A new spherically symmetric general relativistic 
hydrodynamical code. 



Jose V. Romero 1 , Jose M— . Ibanez 2 
Jose M^. Marti 2 ' 3 , and Juan A. Miralles 2 

1 Departamento de Fisica Teorica 

Universidad de Valencia 
46100 Burjassot (Valencia), Spain 
2 Departamento de Astronomia y Astrofisica 
Universidad de Valencia 
46100 Burjassot (Valencia), Spain 
3 Max-Planck-Institut fur Astrophysik 
Karl-Schwarzschild-str., 1, 8046 Garching bei Munchen, Germany 

Abstract 

In this paper we present a full general relativistic one-dimensional hydro-code 
which incorporates a modern high-resolution shock-capturing algorithm, with an ap- 
proximate Riemann solver, for the correct modelling of formation and propagation 
of strong shocks. The efficiency of this code in treating strong shocks is demon- 
strated by some numerical experiments. The interest of this technique in several 
astrophysical scenarios is discussed. 
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I. INTRODUCTION 



Since the pioneering work by May and White (1967) the use of general relativistic and 
spherically symmetric hydro-codes has been restricted, basically, to the field of stellar 
collapse and Supernovae. Currently, there are several astrophysical scenarios for which 
the constraint of spherical symmetry is still a good approximation and where general 
relativistic hydrodynamical processes are involved: i) Gamma-ray bursters: Models for 
explaining gamma-ray bursts compatible with the spatial distribution derived from the 
BATSE experiment from the Compton Observatory are based on relativistic fireballs 
originating from the sudden release of energy in small regions (Meszaros et al. 1993, Piran 
et al. 1993). ii) Spherical accretion onto compact objects: Theoretical studies of spectral 
properties of X-ray radiation produced in atmospheres around an accreting neutron star 
have particular importance given the observational capabilities of present instrumentation 
on board satellites. In particular cases - low magnetic fields, for example - the assumption 
of spherical symmetry is adequate, iii) Stellar collapse: Different topics of current interest 
in the field of stellar collapse and supernovae are, apart from the theoretical problem of 
black hole formation (see, e.g., Baumgarte, Shapiro and Teukolsky 1994), the equation 
of state for dense matter, the role of neutrinos, the influence of convective motions, etc.. 
One of them is the correct modelling of formation and propagation of the shock formed 
after bounce, the so-called prompt phase. This question remains of crucial importance as 
an initial mechanism leading, with the help of other processes involved in the different 
versions of the delayed mechanism - neutrino energy deposition (Bethe and Wilson 1985), 
convective motions (Janka and Miiller 1993, Herant et al. 1994) -, to the final success of 
the explosion. 

The presence of strong shocks is a common feature in the above astrophysical scenarios. 
Shocks are discontinuous solutions of the hydrodynamical equations and are an important 
source of numerical problems and inaccuracies. The correct modelling of strong shocks 
is one of the most delicate issues in - both Newtonian and relativistic - hydrodynamical 
codes. 

May and White's code was built up by using standard finite difference techniques and 
incorporates an artificial viscosity term to damp down spurious numerical oscillations 
around discontinuities. Up to date, codes based on the original formulation of May and 
White and on later versions (e.g., Van Riper 1979) have been used in many supernovae 
calculations (see, e.g., the recent paper by Swesty, Lattimer and Myra 1994 and references 
cited therein). The Lagrangian character of May and White's code together with other 
theoretical considerations concerning the particular coordinate gauge, has prevented its 
extension to multidimensional calculations. 

The first Eulerian code to solve the equations of relativistic hydrodynamics comes from 
the work of Wilson (1972, 1979). The main ideas of Wilson's procedure laid the founda- 
tions of several codes developed in the first half of the eighties which have been applied to 
several astrophysical scenarios: axisymmetric stellar collapse (Piran 1980, Stark and Pi- 
ran 1987, Nakamura et al. 1980, Nakamura 1981, Nakamura and Sato 1982, Evans 1986), 
accretion onto compact objects (Hawley et al. 1984, Petrich et al. 1989) and numerical 
cosmology (Centrella and Wilson 1984). All of these codes make use of a combination 
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of artificial viscosity and upwind techniques in order to obtain numerical solutions of the 
relativistic hydrodynamical equations. On the other hand, the equations are written as 
a set of advection equations, so the terms containing derivatives (in space or time) of 
the pressure are treated as source terms. This procedure breaks down (numerically) the 
conservative character of the relativistic hydrodynamics system of equations (see below). 

Concerning the present work, during the last decade a number of new shock- capturing 
finite difference approximations have been constructed and found to be very useful in the 
numerical simulation of classical (Newtonian) fluid dynamics (see, e.g., LeVeque 1992). In 
addition to conservation form, these schemes are usually constructed to have the following 
properties: a) Stable and sharp discrete shock profiles, b) High accuracy in smooth regions 
of the flow. Schemes with these characteristics are usually known as high-resolution shock- 
capturing schemes (henceforth HRSC). They avoid the use of artificial viscosity terms 
when treating discontinuities and, after extensive experimentation, they appear to be a 
solid alternative to classical methods with artificial viscosity. 

Our spherically symmetric general relativistic hydro-code is the cumulative result of 
our experience in HRSC schemes applied to Newtonian and special-relativistic fluid dy- 
namics (see Ibahez 1993, for a review). In previous papers (Marti et al. 1991, Marquina 
et al. 1992) HRSC schemes have been extended to solve the relativistic hydrodynamics 
system of equations in one spatial dimension. The procedure relied on two main points: 
1) To write the equations of relativistic hydrodynamics as a system of conservation laws 
and identify the suitable vector of unknowns. 2) An approximate Riemann solver built up 
from the spectral decomposition of the Jacobian matrix of the system at the boundaries 
of each numerical cell. Schneider et al. (1993) have also explored the ultrarelativistic 
regime with a different Riemann solver. The multidimensional extension of our special- 
relativistic hydro-code can be found in Font et al. (1994). Similar results have been 
independently derived by Eulderink (1993). Recently, Marti and Muller (1995) have built 
up a relativistic version of the popular Newtonian Piecewise Parabolic Method (PPM, 
Colella and Woodward 1984) and have carried out simulations of relativistic jets with 
Lorentz factors greater than 20 (Marti et al. 1995a). 

Finally, note that in recent years, techniques other than the standard finite difference 
ones have been applied to simulate relativistic flows. Gourgoulhon (1991,1992) made use 
of the accuracy of spectral methods (Gottlieb and Orszag 1977) to detect the zero value of 
the fundamental mode against radial oscillations of a neutron star at the maximum of the 
mass-radius curve for a given equation of state. Although the global error on the solution 
decreases exponentially with the number of degrees of freedom, the handling of shock 
waves by spectral methods remains a major problem due to the Gibbs phenomenon. By 
combining moving grids and shock tracking techniques Bonazzola and Marck (1991) have 
obtained promising results for one-dimensional flows. A different numerical approach is 
provided by smooth particle hydrodynamics (Gingold and Monaghan 1977, Lucy 1977), 
the first application to simulate relativistic flows being made by Mann (1991). These 
results, although encouraging, are still far from those obtained with HRSC methods. 

General relativistic hydrodynamics links gravitational field with geometry, hence, from 
the numerical point of view, new difficulties arise due to the nonlienarities introduced 
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by the geometrical terms. The system of equations is coupled not only by Lorentz-like 
factors, as in special relativity (see Norman and Winkler, 1986 for a discussion on the 
numerical difficulties), but also by the different components of the metric tensor. In the 
present paper we describe the features of a full general relativistic hydrodynamical code 
(some preliminary results were presented in Ibanez et al. 1992) having the following main 
properties: i) The code makes use of an approximate Riemann solver which allows to 
capture shocks in a consistent way. It is a natural extension to general relativistic fluid 
dynamics of modern HRSC schemes, ii) By using splitting techniques the code can be 
easily extended to the multidimensional case once the behaviour of the characteristic fields 
is known. 

The structure of this work is as follows: Section §// is devoted to display the full 
equations of general relativity in a spherically symmetric space-time as well as to the 
theoretical analysis of this system considered as a system of conservation laws. The 
hydro-code is based on a numerical algorithm which is explained in §///. Some numerical 
tests and astrophysical applications are shown in §IV. Finally, a summary of our main 
conclusions is presented in §V. 

II. GENERAL RELATIVISTIC EQUATIONS IN SPHERICAL 
SYMMETRY AS A SYSTEM OF CONSERVATION LAWS 

Let M. be a general spacetime, described by the four dimensional metric tensor g^ u . 
According to the {3 + 1} formalism, the metric is split into the objects a (lapse), (5 l (shift) 
and 7y, keeping the line element in the form: 

ds 2 = -(a 2 - fiip^dt 2 + 2(3 i dx i dt + -f ij dx i dx j (1) 

where Greek (Latin) indices run through all (spatial) coordinates. 

In the case of spherically symmetric space-times the general relativistic equations 
can be given in a simple way which looks like the Newtonian hydrodynamics. To this 
aim the choice of coordinates is crucial. Schwarzschild-type coordinates (Bondi 1964) 
allow a simple extension of the Eulerian Newtonian hydrodynamics to the Einstenian 
one. In terms of slicing of space-time, Schwarzschild-type coordinates are the realization 
of a polar time slicing, a radial gauge (see Gourgoulhon 1991), and the generalization 
of the Schwarzschild coordinates, in terms of which the vacuum and static space-time is 
described. 

The choice of the radial gauge leads to the following expression for the 3-metric 

7jJ =diag(X 2 ,r 2 ,r 2 sin 2 fl) . (2) 

The choice of the polar slicing condition, in spherical symmetry and the radial gauge 
implies a zero-shift vector (f3 % = 0,Vi). In the following, we will use the acronym RGPS 
(radial gauge and polar slicing) for the particular spacetime in which we are interested: 

ds 2 = -a 2 dt 2 + X 2 dr 2 + r 2 (d6 2 + sin 2 6d(f) 2 ) . (3) 
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By analogy with the well-known Schwarzschild solution for vacuum, we can define the 
functions $(r, t) and m(r,t) by 

X(r,t) = (l - ^M^) ' , a (r,t) = exp{$(r,t)} . (4) 

The equations describing the evolution of matter are the expression of the local con- 
servation of baryon number 

V M J" = , (5) 
and the local conservation of energy-momentum 

V M T^ = , (6) 

The current J M and the energy-momentum tensor T^ u are 

J» = pu» , (7) 

Tp, v = phu^u u + pg^ , (8) 

where V M is the covariant derivative, the four-velocity of the fluid, p the rest-mass 
density, p the pressure, and h the specific enthalpy defined by h = 1 + e + p/p, where e 
is the specific internal energy. 

The expression chosen for the energy-momentum tensor is that of a perfect fluid and 
we therefore ignor effects due to heat conduction or viscous interactions. 

Let us introduce the physical velocity, v, defined by v — Xu r /au l . This quantity 
represents the fluid velocity relative to an observer at rest in the coordinate frame. The 
Lorentz-like factor W = av}, satisfies the familiar relation with v, W — (1 — v 2 )^ 1 ^ 2 . 

With the above coordinate conditions and the following set of unknown variables 

D = aXJ t = XpW , 
S = aT tr = phW 2 v , 

r = a 2 T tl — D = phW 2 — p — D , (9) 

the general relativistic equations can be written as a system of conservation laws (with 
sources) : 

where 

u = (D,S,t) , (11) 

is the vector of unknowns which define the state of the system. The fluxes, f , are defined 

as 

f = (Dv, Sv + p, S - Dv) , (12) 
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and the source terms (free of derivatives of hydrodynamic quantities) are 



s(u)= (o,(Sv-r- D)(8aXnrp + aX^-) + aX P ^- + ^-,o] . (13) 
V r 2 r 2 Xr J 

The above system of equations is closed by an equation of state (EOS), that we will 
assume of the form 

P = p(p,t) , (14) 

and the Einstein equations which furnish conditions on the quantities m(r,t) and $(r, t) 
(see, Gourgoulhon 1991): 

QflfT 

_ = 4vrr 2 (r + D) , (15) 
or 



^l = X 2 (^- + A-Kr{p+Sv)) , (10) 



Indeed, by analogy with the static case, we can distinguish (for a given spherical 
surface having an area 4irr 2 ) between the enclosed gravitational mass m defined by fll5D 
and the enclosed baryonic mass vha defined by 

°^ = A^D . (17) 
or 

Equations flI5p and ([T7]) are integrated, at each time step, between r = and r = 
R(t) (the radius of the star) with the following boundary conditions: m(r = 0) = 0, 
m^ir = 0) = 0, m(R(t)) = M (the total gravitational mass), rriA{R{t)) = Ma (the total 
baryonic mass). The gravitational potential $, given by equation ([R]), is defined with the 
exception of an aditive arbitrary constant and is matched at the surface with the exterior 
Schwarzschild's solution, i.e., &(R(t)) = (l/2)ln(l— 2M/R(t)). The integration of $ starts 
with a zero value at the center and by comparing, at the surface, the integrated value of 
$ with the matching condition it is possible to obtain the arbitrary aditive constant. 

It is worth pointing out that the system of equations ( |10D displays a very important 
feature, that is, the conservation of baryonic mass and energy. In effect, the first and 
third equations of system (|T0|) lead (when acting on integral quantities) to the following 
relations: 



dM A a 



-—Dv 



dt X 
dE b a 

~dT ~X 



R(t) 



— (S-D)v \ R[t) 



The quantity E b = M — Ma, on the analogy of the static case, can be considered as the 
binding energy of the star. The symbol |#( t ) means that the quantities on the right hand 
side of these equations are evaluated at the surface. Hence, if the boundary conditions at 
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the surface are zero (i.e., p = e = 0, or v = 0), the conservation of Ma and Ej, (or M) is 
strictly satisfied. 

A numerical algorithm written in conservation form, like the one used in the present 
paper (see below), preserves numerically this important property of the system. Hence, 
the correct choice of the vector of unknowns u, and consequently, source terms s(u) is of 
crucial importance. 

From the conserved quantities C = {D, S, r} we must obtain the set of primitive 
variables p = {p, v,e}, at each time step, by solving an implicit equation in pressure. A 
one-dimensional Newton- Raphson routine suffices to obtain p (see Marti & Miiller 1995, 
for details). In the Newtonian limit, the set of new variables C = {D, S, r} tends to the 
set {p, pv, e}, with e = pe+ (l/2)pv 2 — m/r, i.e., the density, momentum density and total 
energy density, respectively. 

The hyperbolic character of the relativistic hydrodynamic system of equations has 
been the subject of study of many authors (see Anile, 1989 and references therein). 

As we will show in the next section, Godunov-type methods, or modern HRSC schemes, 
incorporate the resolution of local Riemann problems, initial value problems for system 
(|i~0|) with discontinuous data. In order to extend HRSC to our problem is crucial to know 
the spectral decomposition of the Jacobian matrix £>(u) of the system (JTO|) : 

B=^M . (19) 

Following a procedure similar to the one described in Font et al. (1994), we have 
derived the eigenvalues and right-eigenvectors of B(u). 

The eigenvalues, the characteristic speeds associated with material waves and acoustic 
waves, are respectively, 



and 



X = v , (20) 



A ± = ^ , (21) 
1 ± vc q 



where c s is the local sound velocity which satisfies 

hc 2 s = X + (p/p 2 )K , (22) 

with x — (9p/dp) e and k = (dp/de) p . 

These eigenvalues are the natural extension to our generic space-time of the corre- 
sponding characteristic speeds well-known in Minkowski space-time or in Newtonian fluid 
dynamics. Let us point out that they seem like in Minkowski space-time due to our 
definitions of velocity and flux. 

The right-eigenvectors are: 
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where k = k/ p, and r^ ± is the first component of the corresponding eigenvector. 

The above system ([TIT) is strictly hyperbolic since the Jacobian matrix B(u) has real 
and distinct eigenvalues. 

III. OUR MODERN HIGH-RESOLUTION SHOCK-CAPTURING 

SCHEME 

In order to exploit numerically the conservative character of the system fllPP we have 
written it as a hyperbolic system of conservation laws (see Lax 1972, for a mathematical 
analysis). Let us summarize in this Section the main features of our algorithm (see also 
Marti et al. 1991). 

At each time level, the data are the cell averages of the conserved quantities 

u" = ^— f r ' J+h u(r,t n )r 2 dr . (25) 

The data are advanced in time according to a version of the method of lines: 
dt TjAr-j 

where 

f i+ i =f(iif_i, Uj,u j+1 ) , (27) 

is a consistent numerical flux vector, i.e., f(u, u, u) = f(u). The quantity Aj + i is a 
combination of the geometrical factors 



■jtL • (28) 

evaluated at the interface j + h, and s 3 - is the cell average of the source terms calculated 
according to (f25|). 

Once the procedure to evaluate f- + i is known , then system (|26|) can be integrated 
in time by using a suitable ordinary differential equation solver. We have made use of 
a standard predictor-corrector method. The value of the timestep is constrained by the 
Courant condition. 

A reconstruction procedure of the solution at the time level t n from its cell averages 
allows us to define local Riemann problems at each interface j + |. We have used a 
monotonicity preserving linear reconstruction of the primitive variables using the minmod 
function as a 'slope limiter' (Van Leer, 1979). Accordingly, the corresponding values of 
u, i at the interface j + ^are: 

3 + 2 J 2 

u£i= u? + S](r j+h - rj ) , (29) 



S 



where S™ is a s/ope limiter defined by 



u f + i = i ? + i + s iAr r i + i) > (3°) 



/fi« — fi n fi n — ii n \ 
S" = minmod j+1 j , j , (31) 



and the minmod function makes a choice of the slope which is minimum or takes a zero 
when they have different signs: 



minmod(a, b) 



a if | a | < | b | , ab > 
6 if | a | > | b | , ab > 
if ab < 



At each cell interface the i component of the numerical flux is computed as follows 
% = \ ( / (0 «i) + / W « + i) - E I Aa I Atf a ff) , (32) 

\ a=0,± / 

where L and R stand for the left and right states at a given interface j + |, A a and 
(a = 0, ±) are respectively, the eigenvalues and the i i/l -component of the a-right 
eigenvector of the Jacobian matrix 



- V ^ /u-J(u« +i+ u; + ,, 



The quantities Au a , the jumps in the local characteristic variables across each cell 
interface, are obtained from 

uf + i - u^ + , = E A^ a ? a . (34) 

a=0,± 

A a , f a and AcD Q , as functions of u, are evaluated at each interface and, therefore, they 
depend on the particular values u^ +i and u^_i. 

With the above ingredients, our algorithm is conservative, upwind, and due to the 
particular cell-reconstruction, monotone. It is globally second order accurate, although 
this statement is only well-founded for scalar equations, equally spaced grids, and the 
smooth part of the flow. 

Finally, let us comment on one of the features of our code, its modularity. It has been 
constructed to allow for easy substitution of different Riemann solvers and different cell 
reconstructions. It also works in Newtonian and special relativistic hydrodynamics. 

IV. NUMERICAL EXPERIMENTS AND ASTROPHYSICAL 

APPLICATIONS 
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Modern numerical codes must be, and in fact are, checked against well-known an- 
alytical solutions before their exploitation in applications. A long battery of test-beds 
addressed to check hydrodynamical codes exists. The most popular ones are the standard 
shock problems in different symmetries, planar (shock-tube tests), cylindrical or spherical. 
Their corresponding analytical solutions (in Newtonian fluid dynamics) can be found in 
standard textbooks (see,e.g., Courant and Friedrichs 1948). The one-dimensional solu- 
tion of the Riemann problem for relativistic hydrodynamics has been derived recently by 
Marti and Miiller (1994) allowing for the analytical solution of any initial value problem. 
In previous papers we have carried out numerical experiments in planar symmetry with 
the following standard shock-tube problems: Sod's test, blast wave and shock reflection 
tests, in both Newtonian and special-relativistic hydrodynamics (see Marti et al. 1991, 
and Marquina et al. 1992). Here, in Section §/VT, we have concentrated on the relativis- 
tic spherical shock reflection problem, which allows to check special relativistic dynamics 
in a spherically symmetric geometry. 

A second test (see Section %IV.2) is the spherical accretion onto a self-gravitating 
object, including general relativistic effects. This test has the advantage of allowing one 
to check stationary solutions of general relativistic hydrodynamics on a fixed background. 

A third test (see Section §IV.3) underwent by our hydro-code is the detection of 
the zeros of the fundamental mode against radial oscillations of a spherical equilibrium 
configuration, which checks the stability of equilibrium solutions. 

In Section §IVA, we have analyzed Oppenheimer-Snyder collapse. This is one of the 
standard tests for spherically symmetric, fully relativistic codes, which allows one to check 
general relativistic dynamics. 

With the above selection of tests we span a broad range of different features: special 
relativistic effects, geometrical effects, general relativistic effects in a background and fully 
general relativistic dynamics. 

Finally, in Section §IV.5 we have carried out an analysis of the dynamics of the general 
relativistic collapse of compact objects. 

IV. 1 Relativistic spherical shock reflection 

Noh (1987) has exhaustively analized the spherical shock reflection problem in the 
Newtonian dynamics of ideal gases (p = (T — l)pe). Noh's paper is, in fact, devoted 
to the study of intrinsic errors which appear when using artificial viscosity techniques in 
spherically symmetric applications. We have considered the relativistic version of this 
numerical experiment. The relativistic spherical shock reflection is a severe test due to 
the difficulties connected with the geometry and, in its relativistic version, with strong 
nonlinearities induced by the Lorentz factor in the ultrarelativistic regime. 

Analytical solutions for the relativistic shock reflection with planar, cylindrical and 
spherical symmetry can be found in Marti et al. (1995b). The initial data consists of 
an inflowing cold (i.e., e = 0) gas with coordinate velocity v and corresponding Lorentz 
factor W. In the spherical case, the flow converges towards the center and its reflection 
causes the compression and heating of the gas as it converts its momentum into internal 
energy. A shock is formed which starts to propagate through the inflowing gas. Behind 
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the shock the gas is at rest and, according to the conservation of energy accross the shock, 
has a specific internal energy given by 



e+ = W-l. (35) 

The jump in density at the shock is defined by Ap = p + / p~ , where p + and p~ are the 
postshock and preshock values, respectively. For a relativistic strong shock, like the one 
developed herein, this jump satisfies 

YW + 1 

Ap = T3 ^, (36) 



with a shock velocity given by 



(r- \)w\v\ 



- w+ i ■ (37) 

In the nonshocked part of the flow (i.e., r <E]v s t,oo[), the proper rest-mass density 
distribution is given by 

p = (l + jj p (38) 

where po is the rest-mass density at infinity. 

It is worthwhile to point out one relevant difference between classical and relativistic 
shocks. In the relativistic case, the jump in density is unbounded, whereas in the present 
test, the ultrarelativistic regime leads to the following asymptotic relations (p = 1), 

Ap -> — !— W , 
' T- 1 



P + 



r 



i) w • 



r 

(r-i) . (39) 



The initial data for this problem are defined in a unit sphere (0 < r < 1) and for an 
ideal gas with T = 4/3; v(r, 0) = —vo(vo > 0), p(r, 0) = po = 1, e(r,0) = eo = 0. For 
numerical reasons, the initial value of the specific internal energy of the inflow gas was set 
to a small nonzero value e = lO -6 !^- The boundary conditions are v(l,t) = —v and 
the exact solution of p at r = 1 for all time. 

We switched off the gravitational terms in our code and transformed it into a purely 
special relativistic hydrodynamical code, and have used a grid with 200 equidistant zones. 
The results displayed in Table 1 and Fig.l, correspond to a particular instant of the 
evolution after 1000 timesteps. In order to emphasize the relativistic effects, we have 
done the calculation for a large sample of different initial inflow velocities Vq. Several 
conclusions can be established from the data contained in Table 1: i) The postshock 
density increases with the initial inflow velocity and is unbounded. This is a typical 
relativistic effect, ii) The ratio p + /Wo tends to the asymptotic value of 64 when t> — > 
1, which is consistent with the above asymptotic relations, iii) The maximum of the 
relative errors for the postshock density, e™ ax , and the mean relative error of the same 
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quantity, e r , are independent of the initial Lorentz factor and converge to 14% and 2%, 
respectively. This feature of relative errors has been pointed out by Marti and Miiller 
(1995) in discussions of planar shock wall test experiment with their relativistic version 
of PPM. In the mean relative error we have not considered the zone next to the center, 
which always dominates the maximum error, due to the well-known effect of numerical 
overheating. Figure 1 shows our numerical results compared with the analytical ones. 
The overheating phenomena is less severe by far than the one found by Noh (1987) in 
his analysis of the corresponding Newtonian problem. Noh (1987) reported errors of the 
order of 1000% in some of the experiments (in Newtonian hydrodynamics) with artificial 
viscosity techniques. A maximum relative error of 14% in the postshock density (at the 
center) is comparable to the one obtained with the Newtonian PPM reported by Noh 
(1987, see Fig. 24). 

IV. 2 Spherical accretion onto a black hole 

We have tried to reproduce some of the stationary solutions of the spherical accretion 
onto a black hole in two cases: i) dust accreting onto a Schwarzschild black hole (Hawley 
et al. 1984) and ii) an ideal gas accreting onto a Schwarzschild black hole (Michel, 1972). 

In these tests, the gravitational field is kept fixed. The initial conditions are those 
of a vacuum, i.e., density is zero everywhere except at the outer boundary where a gas 
is being continously injected with a velocity and density given by the exact solution. 
We have run our code using two different grids of 50 and 100 points, which span the 
interval 1.05 < r/2M < 10.0. Outflow boundary conditions have been taken at the inner 
boundary r = 2.1M. 

The analytical solution of a spherical geodesic (presureless) flow accreting onto a black 
hole can be found in Hawley et al. (1984). Table 2 summarizes relative errors in the three 
fundamental variables when compared with the analytical solution at a time of 180 M 
(when the stationary solution has been reached). Discrepancies with the exact solution 
amount to less than 2% for the maximum of the relative errors and less than 1% for 
the mean relative error. As it should happen these errors decrease under grid refinement 
(compare first and second rows in Table 2). 

The stationary solution of an ideal gas accreting spherically onto a Schwarzschild black 
hole was derived by Michel (1972). A critical point exists, as in the Newtonian description, 
which can create numerical difficulties. We have considered only solutions of accreting 
gas having a critical point far outside of our computational domain. Table 2 summarizes 
the relative errors in the three fundamental variables when compared with the analytical 
solution at a time of 720 M, when the stationary solution has been reached. Discrepancies 
with the exact solution amount to less than 3% for the maximum of the relative errors 
(less than 1% for the mean relative error). As before, these errors decrease when the grid 
is refined. Figures 2 and 3 show relativistic rest mass density and the velocity profiles. 
Three snapshots of quantity D are plotted in Fig. 2 and show that the numerical solution 
converges to the stationary analytical one. The velocity converges so fast that it's not 
possible to distinguish the same snapshots in Fig. 3. 
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As can be deduced from the behaviour of the mean errors under grid refinement (see 
Table 2), our code is second order accurate for continuous solutions. 



IV. 3 Detection of the zero value for the fundamental mode against radial 

oscillations of a compact object 

According to the static stability criterion (Harrison et al. 1965) it is possible to 
establish a correlation between the critical points of a curve (in a "gravitational mass 
versus radius" M — R diagram) built up from equilibrium configurations obeying a given 
cold EOS and the onset of instability of the zero- frequency mode against radial oscillations. 
This criterion avoids solving the Sturm-Liouville equation governing radial perturbations 
of an equilibrium model (see Shapiro and Teukolsky 1983, for a general discussion, or 
Marti et al. 1988, for details on the numerical solution) A hydro-code should allow for 
the dynamical study of perturbations over an equilibrium configuration which can be 
compared with predictions of the static stability criterion. The dynamical detection of 
a zero in the fundamental mode has been used as a test-bed of hydro-codes by several 
authors (see, e.g., Ibanez 1984 and Gourgoulhon 1991). The most accurate results known 
to the authors are those obtained by Gourgoulhon (1991) using pseudospectral techniques. 

In our case we have focussed on a particular critical point, the maximum of the M — R 
curves corresponding to two kind of compact objects (white dwarfs and neutron stars) 
which obey suitable EOS (see below) and such that their stability properties have been 
exhaustively studied (see references below). That maximum separates the curve in two 
branches, being the stable (unstable) one the corresponding to models having central 
denstities lower (higher) than the critical one. 

The calculations of this Section as well as those of Section §IV.5 have been carried out 
with a numerical grid (which is Eulerian, i.e., fixed) built up in such a way that the radius 
of the initial model is partitioned into N zones, distributed in geometric progression in 
order to have finer resolution near the center. In our simulations we have taken N=330 
and the surface radius R(t) is given by the condition p(R(t)) = p(R(t = 0)) ~ 10 5 g/cm 3 . 
Hence, the data relative to the initial stellar model are given in the entire grid. As time 
goes on and, according to the dynamics of each problem, the number of cells covering the 
star decreases. During the contraction phase of our numerical experiments in this Section 
(§/V.3), or during the infall epoch of the applications in Section §IV.5, the surface is 
moved in accordance with the above prescription (p(R(t) = p{R{t = 0))) and we eliminate 
those cells remaining outside R(t) (in practice, we impose over these cells the vacuum 
conditions). The important point is that in none of these applications the number of 
numerical cells used for describing the star are less than 200 (for N=330) keeping the 
numerical resolution in a good level for an accurate calculation. An analogous procedure 
is established for the expanding phase. Other prescriptions for the surface condition such 
as, e.g., a combination of moving grids with the integration of the velocity at the surface 
would be interesting to be studied. 
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IV. 3.1 Detection of a zero in the fundamental mode against radial oscillations of a 
white dwarf 

We have considered a sample of initial models which are white dwarf-like configura- 
tions, obtained by solving the stellar structure equations with the equation of state of an 
ideal gas of electrons including Coulombian corrections due to the ions (Ibanez, 1984) and 
with a homogeneous chemical composition (Carbon). The maximum gravitational mass 
model has a value of 1.3862M . 

We have found a change in the stability behaviour at some point between 1.3847M 
which is a stable model and 1.3853M Q which is an unstable one and collapses in a time 
which is an order of magnitude higher than the characteristic dynamical one. The cor- 
responding values of the initial central densities are 1.3 and 1.5, respectively, in units of 
10 10 gcm~ 3 . 

Figures 4 and 5 show the velocity field as a function of time and radius for, respectively, 
the stable and unstable models. We point out the different temporal scales involved and 
the different interval of values spanned by the velocity in each case. 

We have also run simulations with a grid of 180 points and taken initial models of 
constant density. These models are not equilibrium models and, consequently, they start 
to collapse until some new equilibrium configuration is reached. We find a change in the 
stability of the models at some point between 1.37M and 1.38M Q which is unstable. 

IV. 3. 2 Detection of a zero in the fundamental mode against radial oscillations of a 
neutron star 

We have considered a family of neutron stars obeying the EOS derived by Diaz-Alonso 
(1985) in a field theoretical model approach to neutron matter at zero and finite temper- 
ature (Diaz-Alonso et al 1989). 

The model describes a many-body system of mutually interacting nucleons and mesons: 
a fictitious a particle, pions, p and uj mesons. The interaction is given by a relativistic 
Lagrangian containing the free Dirac Lagrangian for the nucleons, the Lagrangians for 
the a particle, pion and uj and p mesons. The interaction is described by two pieces: 
one containing the meson-nucleon interaction in the form of Yukawa-like couplings, and 
one which describes the meson-meson interaction. The attractive part of the nuclear 
interaction is provided by the a and pion exchange. The uj and p mesons give rise to a 
repulsive short-range interaction, which is charge- dependent in the latter case, in such a 
way that the EOS behaviour is different for nuclear symmetric matter as for pure neutron 
matter. Details on the solution of the model can be found in the references above. 

We will focuss on the EOS II in Diaz-Alonso (1985) which was derived fitting the free 
parameters in the Lagrangian in such a way that, for symmetric nuclear matter, the model 
saturates at a density p n = 2.837 x 10 14 gcm~ 3 . The corresponding binding energy 
is —15.68 MeV, the symmetry energy 33.17 MeV, and the nuclear incompressibility is 
K = 225 MeV. 

The equilibrium configurations obeying this EOS at zero temperature have been taken 
as initial models in our hydro-code. Their macroscopic properties were exhaustively stud- 
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ied in Diaz-Alonso and Ibanez (1985) for the cold version of EOS II. The maximum gravita- 
tional mass is 1.94556M Q which corresponds a central energy density of 2.48 x 10 15 gcm~ 3 . 
The reader interested on the properties of stability against radial oscillations and the 
equilibrium features of slowly rotating stars obeying the hot version of that EOS can 
address, respectively, to Marti et al. (1988) and Romero et al. (1992). 

We have generated a perturbation of the equilibrium models by imposing an initial 
velocity profile according to the law v (r) = —v (r/R). This procedure has the advantage 
of allowing a clear distinction between the dynamics of the models, erasing round-off 
errors induced in the numerical construction of the initial model (for the hydro-code) 
from the model generated by solving the structure equations. As before, the value of R(t) 
is obtained from the condition: p(R(t)) = p(R(t = 0)). 

Taking a value for v Q of 10~ 3 we have found a change in the stability behaviour between 
the model of 1.94532M Q (central density: 2.55 x 10 15 gcm~ 3 ), which is stable, and the 
model 1.94518M (central density: 2.57 x 10 15 gcm -3 ), which is unstable. Unlike the 
stable configuration -which undergoes a oscillating motion- the unstable one collapses 
in a time of the order of the characteristic dynamical time. The behaviour of the rest- 
mass density, as a function of time, can be seen in Figs. 6 and 7 for, respectively, the 
stable and unstable models. In these figures we have selected a sample of mass shells 
by making a equipartition of the total gravitational mass in eleven shells of the same 
width and plotted the ten inner ones. For the values of the central density involved in 
this study, the dynamical characteristic time is of about < 0.8 msec, which is interesting 
in view of the different temporal scales which govern the dynamics of both models. The 
stable configuration (Fig. 6) displays a plateau profile during more than 5 times its 
dynamical characteristic time. The unstable configuration (Fig. 7) starts to collapse in a 
characteristic time less than two times its dynamical characteristic time. 



IV. 4 Oppenheimer-Snyder collapse 

The gravitational collapse (into a black hole) of a homogeneous spherical dust cloud 
(p=0) has been exhaustively studied since the original paper of Oppenheimer and Snyder 
(1939). Depending on the particular time slicing and coordinate gauge used simple ana- 
lytical expressions are available for both the metric coefficients and the matter variables 
(Petrich et al., 1985, 1986). The solution to this problem in the RGPS (as in the present 
work) can be found in Petrich et al. (1986) and Gourgoulhon (1992,1993). We have 
summarized these previous theoretical works in an appendix. Since during the last epoch 
of Oppenheimer-Snyder collapse, the variables involved (both the geometrical and hydro- 
dynamical) reach extreme values, it is generally considered a good test-bed calculation 
for fully general-relativistic time-dependent numerical codes (Petrich et al., 1985, 1986; 
Gourgoulhon, 1992,1993; Baumgarte et al., 1994). Hence, we have used it as a check of 
our code in the case of strong gravitational fields, as well as ultrarelativistic velocities. 
Although no shocks appear during the evolution of the presureless ball, very steep spatial 
and temporal gradients develop at the last stages of its evolution. 
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Results of our simulations have been plotted together with the exact solution in Figs. 
8-10, as a function of radial coordinate, normalized to the Schwarzschild radius of the 
initial configuration, for a sample of values of the temporal coordinate labelled from 
to 6 (see the corresponding captions). Figure 8 shows the rest-mass density, normalized 
to its initial value. Due to the gauge used, the existence of a positive spatial gradient of 
rest-mass density is noticeable. Figure 9 displays the behaviour of the physical velocity. 
As expected, the velocity of the surface tends to unity as the collapse proceeds to its late 
stage, having a minimum value of about v = —0.92 (plotted in Fig. 9). The geometrical 
quantity a is plotted in Fig. 10. At the late epoch we have succeeded in reaching a value 
as low as 1.3232 x 10" 10 . 

In this application the position of the surface R(t) is defined by the analytical solution. 
Note that due to our procedure for moving the surface the value corresponding to this 
point (the last numerical cell) suffers of an error. This error is associated to the fact 
that our grid is Eulerian making difficult to precise the exact position of the surface. A 
comoving grid would allow to solve the surface with more accuracy and would result in 
spatial gradients steeper than those displayed in Figs. 8-10. 

Figures 11 show in a space-time diagram the trajectories of the different mass shells 
in terms of the proper time (Fig. 11a) and the coordinate time (Fig. lib). The effect of 
freezing in the evolution of the system is noticeable from Fig. 11a, by comparison with 
Fig. lib. Our calculation has evolved to an epoch later than that reached by Gourgoulhon 
(1992) with the pseudospectral technique. The results displayed in Figs. 8-11 confirm the 
powerful capabilities of our numerical techniques for treating steep spatial and temporal 
gradients. 

IV. 5 An application to the stellar core collapse 

A very simple way of modelling the essential features of stellar core collapse in massive 
stars is to incorporate a simple equation of state into a hydro-code, like that of an ideal gas, 
but taking; ibatic exponent which can depend on density according to a particular 

prescription. The EOS we have used is a T-law such that T varies with density according 
to: 

r = T min + 7](\ogp - logp 6 ) , (40) 

with: i] = if p < pb and i] > otherwise (Van Riper, 1979). 

Two set of values for the parameters T min , r\ and pb have been considered: {1.33, 1, 
2.5xl0 14 gcm~ 3 } (model A) and {1.33, 5, 2.5xl0 15 gem" 3 } (model B). Model A exhibits 
standard values of the parameters, i.e., the effective adiabatic exponent of infalling mate- 
rial and the value of nuclear density (for symmetric nuclear matter at zero temperature) at 
the saturation point. Model B is rather exotic due to the particular values for the bounce 
density and stiffness rj = 5; however, we have considered this model in order to check the 
ability of our hydro-code in solving numerically flows which develop strong shocks in very 
strong gravitational fields. 

The initial model in the present application is a white dwarf having a gravitational 
mass of 1.3862M Q corresponding to the maximum mass cited in the section before. 
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We have run models A and B with a grid of 330 points distributed - as in the previous 
application - in a geometric progression (in the following, figures are numbered with a 
number followed by a or b corresponding, respectively, to cases A and B). The solution 
converges quickly by taking different grids. No relevant differences were found in the 
numerical results when the number of points is increased from 330 points to 630 grid 
points. Poor grids lead to results which differ from the converged ones, but the essential 
features of the shock remain sharply solved, in typically two numerical cells. This property 
is one of the most relevant of these numerical techniques. As in Section §IV.3 the value 
of R(t) is obtained from the condition: p(R(t)) = p(R(t = 0)). 

Table 3 shows the main features of both models during the infall epoch. Important 
numbers are those corresponding to the maximum of velocity, that is, 0.41 and 0.62, 
respectively, for cases A and B (see also Figs. 12). Also, it is interesting to point out the 
particular values of the geometric factor a 2 at the maximum compression and to compare 
them with those corresponding at the surface of a typical neutron star ps 0.75. 

The kinetic energy of the material ejected has been calculated at different times, those 
corresponding to some fixed values for the position of the shock in our Eulerian frame. 
Table 4 displays the main features of both models during the prompt phase. Case B is 
globally more energetic than case A. The ulterior evolution after the shock has arrived at 
the surface has not been followed up in the present work. 

Note the behaviour of the velocity field as it can be seen in Figs. 12. The shock is 
sharply solved in one or two zones and is free of spurious oscillations. In these figures, 
labels stand for the temporal sequence of each curve (see the corresponding table caption). 
From Figs. 12 we can see that the radius of the inner core, at the time of maximum 
compression and for which the infall velocity is maximum (see Table 3) and the shock has 
been formed, is 12.6 km. At the corresponding time in case B, the size is only 6.3 km. 
This strong difference together with other elements such as the particular values of the 
velocity and the low values of the geometric factor a 2 (see Table 3) are the signatures of 
the fact that general relativistic effects are more important in case B than in case A. 

Rest-mass density, as a function of radial coordinate, for several times is plotted in 
Figs. 13. The propagation of the shock can be followed in this figure. Qualitatively, case 
B displays the same behaviour. The jump in density is of about one order of magnitude. 
The maximum values of the central density differ in both cases by a factor of five (see 
Table 3). 

Figures 14 show several snapshots of the internal specific energy (in units of c 2 ) as 
a function of radial coordinate. The shock is sharply solved and its propagation is well 
defined there. Note the huge values of the internal energy reached at the center (curve 
labelled by 3 in Figs. 14) of ~ 0.18 (case A), being ps 0.42 in case B. 

The conservative features of our hydro-code, consistent with the conservation laws of 
baryonic mass and gravitational mass (or binding energy), are displayed in Figs. 15 and 
16. Figures 15 show several snapshots of the gravitational mass, as a function of radial 
coordinate. The binding energy has been plotted, as a function of radial coordinate, in 
Figs. 16. In both figures, curves labelled by 1 correspond to the initial model, and the 
ones labelled by 3 to the time of maximum compression. A glance at these figures gives 
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confidence in the conservative features of the code, since the relative errors at the surface 
are consistent with the accuracy of our algorithm. Let us focus on those curves labelled 
by 3: Case A generates an inner core (the seed of the protoneutron star) which has a 
mass and radius greater than those corresponding to case B. The size in mass of the 
inner core at the time of maximum compression is ~ 1.15M Q and ~ l.OM for cases A 
and B, respectively. The total binding energy of the initial model is, in units of M Q c 2 , 
—3.1 x 10~ 4 . Conservation of binding energy is preserved by our code. 

We have built up spacetime diagrams in order to simplify the understanding of the 
evolution and, eventually, to compare with other calculations. Figs. 17 show the areal 
radii which enclose a sample of mass shells (see captions of these figures) in terms of the 
time coordinate. Unlike the Oppenheimer-Snyder case (Fig. 11a), in Figs. 17 is very 
remarkable the existence of an absolute minimum of the radius enclosing a given mass, as 
well as the propagation of the shock. From these figures we can distinguish between the 
inner core which reaches hydrostatic equilibrium and the outer part which is ejected with 
a kinetic energy greater than 4 x 10 51 erg (7 x 10 51 erg) in case A (B) when the shock is 
at 560 km (see table 4). 

Finally, to emphasize the importance of general relativistic effects we plot (in Figs. 
18) the geometrical quantity a 2 as a function of radial coordinate and at several times of 
its evolution. As can be seen from Figs. 18, a 2 is a continuous function throughout the 
shock. The curve labelled 2 corresponds to the time at which the absolute minimum is 
reached at the center, being 0.49 and 0.14 for cases A and B, respectively. 

VII. CONCLUSIONS 

In this work we have described a full general relativistic one-dimensional hydrody- 
namical code which incorporates a modern high-resolution shock-capturing algorithm for 
the correct modelling of formation and propagation of strong shocks. Strong shocks are 
sharply solved. Our algorithm is conservative, monotone and upwind. It makes use of a 
linearized Riemann solver. The present version of our hydro-code has the fundamental 
property of conservation of those quantities (such as the baryonic mass and the total 
energy) whose evolution is described by continuity-like equations. 

We have carried out several numerical tests and applications of our code. They have 
been selected in order to check a large range of different properties: special relativistic 
effects, geometrical effects, general relativistic effects in a background or fully general rela- 
tivistic dynamics. In particular, the spherical shock reflection problem has been solved in 
the ultrarelativistic regime most successfully. We have reproduced some of the stationary 
solutions describing a flow evolving in a given background. We have compared dynamical 
study of perturbations with the static stability criterion of equilibrium configurations. 
Oppenheimer-Snyder collapse has given the opportunity to check our code with the ana- 
lytical solution of a fully dynamical spacetime. Finally, we have made a simple application 
to the dynamics of the collapse of compact objects using a simple microphysics. Case B of 
this last application displays strong shocks evolving in presence of gravitational fields so 
huge that the coupling introduced by the geometrical quantities makes the system to be 
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solved a highly nonlinear system. Hence, the severity of some of the tests that our code 
has overcome lead us to be confident in the quality of the results in future applications. 

Let us give some examples of the astrophysical applications that we are envisaging. 
First, we are interested in carrying out simulations of collapsing stellar cores with a 
realistic initial model and an updated microphysics. The influence of the gravitational field 
in the formation and propagation of relativistic fireballs - considered as good candidates 
for modelling 7-ray bursts - will also be studied with the hydro-code analyzed in present 
paper. 
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Appendix A. Oppenheimer-Snyder collapse in RGPS coordinates 

The gravitational collapse of a homogeneous and presureless ball of dust, initially at 
rest (Oppenheimer and Snyder, 1939), has been proposed as a test-bed for fully general- 
relativistic codes. Let us summarize in this appendix the analytical results of this problem 
when they are expressed in terms of the RGPS coordinates (Petrich et al. 1986, Gour- 
goulhon 1993). 

According to Petrich et al. (1986), the solution consists of finding coordinate trans- 
formations from the canonical forms of the interior Friedmann -closed- and the exterior 
Schwarzschild metrics. 

The closed Friedmann metric is 



where r is the proper time of the fluid and x varies into the interval 0<x<Xs,X = 
and Xs being, respectively, the values of x & t the center and at the surface of the star. Xs 
can be determined from the gravitational mass of the cloud M and its initial radius R(0): 
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ds 2 = -dr 2 + a(r) 2 (d X 2 + sin 2 X (d6 2 + sin 2 6d(f) 2 )) , 



(41) 




(42) 



In terms of the conformal time 77, defined by the relation 



dr] 1 



(43) 



dr a(r) 
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the metric (IT) can be written 

ds 2 = a{ V ) 2 (-drj 2 + d X 2 + sin 2 X {d6 2 + sin 2 6d<f) 2 )) . (44) 
The exterior metric is given by the well-known Schwarzschild solution for vacuum 



ds z 



2M\ „ 
1 )dt 2 



2M 



dr 2 + r 2 (d9 2 + sin 2 



(45) 



The RGPS form of the interior metric is given by Eq. (||]). 

The solution of the Oppenheimer-Snyder in the RGPS spacetime can be derived in, 



basically, two steps: i) First, the interior Friedmann metric fl44|) has to be transformed into 
RGPS coordinates (|3]). ii) Second, interior and exterior solutions have to be appropriately 
matched at the surface of the star. 

The explicit expressions we were looking for are: 

1) Coordinate transformations: 



t = 2M 



if 1 A /tan(?? s /2) — tanx 



tanx 

where rj s = r) s (r] c ) = — 2arccos 



2 sin Xs 

' cos(ri c (t)/2) \ 



(46) 



+ tan x s , 

and rj c = r) c (t) being the value of r\ at 



cos 1 / 2 Xs ) 

X = on the hypersurface S t (t=constant). The above equation allows to obtain -by 
inverting it- the function i] c (t). 

The radius at the surface satisfies: 



\ COSXs 



The function %(t, r) is implicitly defined by: 



r = R(0) M 
2) Geometric quantities: 



cos 



'\Vc(t)/2) \ 



smx 



cosx J smxs 



a(r,t) 



cos y — cos 2 (r7 c /2) „ . 



(cos 3 x — cos 2 '{j]c /2)) 1 ' 12 



2M\ 1/2 



if r > R(t) 



(47) 



(4£ 



X(r,t) 



cosx — cos 2 (r/ c /2) 
cos 3 x ~ c °s 2 ('7c/2) / 

2M\" 1/2 



1/2 



if < r < i?(t) 



if r > 
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where 

( \ = cos 3 Xs ~ cos 2 (?? c /2) 
(cos Xs - cos 2 (r] c /2)y/ 2 



3) Matter variables: 



v(r t) - cos(77 c /2)tan X 

(cosx-cos 2 (r/ c /2))V2 ^ 

„( r t ) = 3 sin6 I cos * ^ 3 f5n 

^ ' ; 32ttM 2 l.cosx-cos 2 ^)^ 1 J 
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Table 1 

Relativistic spherical shock reflection 



U \J 


Wn 


P + 

r 


_max 




0.1 


1.005 


342 


0.088 


0.016 


0.9 


2.3 


343 


0.090 


0.015 


0.99 


7.1 


614 


0.11 


0.018 


0.999 


22 


1580 


0.13 


0.020 


0.9999 


71 


4671 


0.14 


0.021 


0.99999 


224 


14455 


0.14 


0.021 


0.999999 


707 


45399 


0.14 


0.022 


0.9999999 


2236 


143252 


0.14 


0.022 



Initial inflow velocity is —vq (in units of the speed of light) corresponding to an initial 
Lorentz factor W . Maximum and mean relative errors of the postshock density, e™ ax 

and e r (after 1000 timesteps). 



Table 2 

Spherical accretion onto a black hole 







D 


s 


T 


EOS 


Grid 






^max 


Dust 
(p=0) 


50 
100 


0.020 0.008 
0.006 0.002 


0.020 0.009 
0.006 0.002 


0.022 0.010 
0.007 0.003 


Ideal 

gas 


50 
100 


0.021 0.009 
0.006 0.002 


0.030 0.010 
0.009 0.003 


0.033 0.012 
0.010 0.003 



Maximum and mean relative errors of the numerical solution, e max and e (at a time 720 

M). EOS: equation of state. 
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Table 3 



Features of the stellar 
collapse: infall epoch 



Model 


Pl4 V max 0t min 


A 


7.28 0.41 0.49 


B 


34.6 0.62 0.14 



p™ x is the maximum of central density at the infall epoch (in units of 10 14 gem -3 ). 
— v max is the maximum of velocity at the infall epoch (in units of the speed of light). 
o^in is the minimum of the purely temporal component of the metric tensor at the infall 

epoch. 



Table 4 



Features of the stellar 
collapse: prompt phase 



Model 


KEi KE 2 KE 3 


A 


1.26 3.18 4.53 


B 


1.18 4.93 7.11 



KEi kinetic energy of the material reaching escape velocities, when the shock is at the 
position Ti (=100, 240, 560 km. Respectively). Energies are given in units of 10 51 erg. 
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Figure captions 

Figure 1.- Relativistic spherical shock reflection. We use a continuous line for the 
exact solution and different symbols like a plus, a diamond and a triangle for velocity (in 
units of the speed of light), density and specific internal energy, respectively. 

Figure 2.- Spherical accretion of dust onto a Schwarzschild black hole. Snapshots of 
the relativistic mass density, in logarithmic scale and geometrized units, as a function of 
the radial coordinate (in units of the mass of the black hole) 

Figure 3.- Spherical accretion of dust onto a Schwarzschild black hole. Velocity versus 
the radial coordinate (in units of the mass of the black hole) at a particular time. 

Figure 4.- Velocity profile as a function of radial coordinate and time for a stable 
white-dwarf. 

Figure 5.- Velocity profile as a function of radial coordinate and time for an unstable 
white-dwarf. 

Figure 6.- Rest-mass density as a function of time for a stable model of neutron 
star having M = 1.94532M and a central density of 2.55 x 10 15 g/cm 3 . Each curve 
corresponds to the following mass shells: rrij = (M/ll) x j (j = 1, ...10) 

Figure 7.- Rest-mass density as a function of time for an unstable model of neutron 
star having M = 1.94518M and a central density of 2.57 x 10 15 g/cm 3 . Each curve 
corresponds to the following mass shells: rrij = (M/ll) x j (j = 1, ...10) 

Figure 8.- Snapshots of the rest-mass density (normalized to the initial value) as a 
function of the radial coordinate (normalized to 2M) for the Oppenheimer-Snyder collapse. 
Labels stand, respectively, for the following values of time (in msec): 0, 0.155, 0.189, 
0.206, 0.218, 0.228 and 0.430 

Figure 9.- Snapshots of the velocity as a function of the radial coordinate (normalized 
to 2M) for the Oppenheimer-Snyder collapse. Labels stand, respectively, for the following 
values of time (in msec): 0, 0.155, 0.189, 0.206, 0.218, 0.228 and 0.430 

Figure 10.- Snapshots of the lapse as a function of the radial coordinate (normalized 
to 2M) for the Oppenheimer-Snyder collapse. Labels stand, respectively, for the following 
values of time (in msec): 0, 0.155, 0.189, 0.206, 0.218, 0.228 and 0.430 

Figure 11a.- Spacetime diagram for the Oppenheimer-Snyder collapse. Trajectories 
of a sample of mass shells in a proper time (in msec.) versus radial coordinate (normalized 
to 2M). Each curve corresponds to the following mass shells: rrij = (M/ll) xj (j = 1, ...10) 

Figure lib.- Spacetime diagram for the Oppenheimer-Snyder collapse. Trajecto- 
ries of a sample of mass shells in a coordinate time (in msec.) versus radial coordi- 
nate (normalized to 2M) diagram. Each curve corresponds to the following mass shells: 
rrij = (M/ll) xj(j = 1,... 10) 

Figure 12a.- Snapshots of the velocity profile as a function of the radial coordinate 
(model A). Labels stand , respectively, for the following values of time (in msec): 80.49, 
80.74, 80.91, 81.04, 81.20, 81.39, 81.76, 81.94, 82.86, 84.25 and 86.12 

Figure 12b.- Snapshots of the velocity profile as a function of the radial coordinate 
(model B). Labels stand , respectively, for the following values of time (in msec): 80.68, 
80.79, 80.93, 81.06, 81.25, 81.33, 81.76, 82.47, 83.22, 84.73 and 86.51 

Figure 13a.- Snapshots of the rest-mass density profile as a function of the radial 
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coordinate (model A). Labels stand , respectively, for the following values of time (in 
msec): 80.49, 80.91, 81.20, 81.76, 82.86 and 86.12 

Figure 13b.- Snapshots of the rest-mass density profile as a function of the radial 
coordinate (model B). Labels stand , respectively, for the following values of time (in 
msec): 80.68, 80.93, 81.25, 81.76, 83.22 and 86.51 

Figure 14a.- Snapshots of the specific internal energy profile as a function of the 
radial coordinate (model A). Labels stand , respectively, for the following values of time 
(in msec): 80.49, 80.74, 80.91, 81.04, 81.20, 81.39, 81.76, 81.94, 82.86, 84.25, and 86.12 

Figure 14b.- Snapshots of the specific internal energy profile as a function of the 
radial coordinate (model B). Labels stand , respectively, for the following values of time 
(in msec): 80.68, 80.79, 80.93, 81.06, 81.25, 81.33, 81.76, 82.47, 83.22, 84.73 and 86.51 

Figure 15a.- Snapshots of the gravitational mass profile as a function of the radial 
coordinate (model A). Labels stand , respectively, for the following values of time (in 
msec): 0, 80.49, 80.91, 81.20, 81.76, 82.86 and 86.12 

Figure 15b.- Snapshots of the gravitational mass profile as a function of the radial 
coordinate (model B). Labels stand , respectively, for the following values of time (in 
msec): 0, 80.68, 80.93, 81.25, 81.76, 83.22 and 86.51 

Figure 16a.- Snapshots of the binding energy profile as a function of the radial 
coordinate (model A). Labels stand , respectively, for the following values of time (in 
msec): 0, 80.49, 80.91, 81.20, 81.76, 82.86 and 86.12 

Figure 16b.- Snapshots of the binding energy profile as a function of the radial 
coordinate (model B). Labels stand , respectively, for the following values of time (in 
msec): 0, 80.68, 80.93, 81.25, 81.76, 83.22 and 86.51 

Figure 17a.- Spacetime diagram for model A. Trajectories of a sample of mass shells 
in a proper time (in msec.) versus radial coordinate (in km. and logarithmic scale). Each 
curve corresponds to the following mass shells: rrij = (M/ll) x j (j = 1, ...10) 

Figure 17b.- Spacetime diagram for model B. Trajectories of a sample of mass shells 
in a proper time (in msec.) versus radial coordinate (in km. and logarithmic scale). Each 
curve corresponds to the following mass shells: rrij = (M/ll) x j (j = 1, ...10) 

Figure 18a.- Snapshots of the geometrical quantity a 2 as a function of the radial 
coordinate (model A). Labels stand , respectively, for the following values of time (in 
msec): 80.49, 80.91, 81.76, and 86.12 

Figure 18b.- Snapshots of the geometrical quantity a 2 as a function of the radial 
coordinate (model B). Labels stand , respectively, for the following values of time (in 
msec): 80.68, 80.93, 81.76 and 86.51 
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IDEAL GAS ACCRETION 
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